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ABSTRACT 

Comptonization is the process in which photon spectrum changes due to multiple Compton 
scatterings in the electronic plasma. It plays an important role in the spectral formation of 
astrophysical X-ray and gamma-ray sources. There are several intrinsic limitations for the 
analytical method in dealing with the Comptonization problem and Monte Carlo simulation 
is one of the few alternatives. We describe an efficient Monte Carlo method that can solve 
the Comptonization problem in a fully relativistic way. We expanded the method so that it is 
capable of simulating Comptonization in the media where electron density and temperature varies 
discontinuously from one region to the other and in the isothermal media where density varies 
continuously along photon paths. The algorithms are presented in detail to facilitate computer 
code implementation. We also present a few examples of its application to the astrophysical 
research. 



1. Introduction 

Comptonization - the process where photon spectrum changes due to multiple Compton 
scatterings in the electronic plasma - is one of the most important processes in the spectral 
generation of X-ray binaries, active galactic nuclei and other X-ray and gamma-ray sources. The 
analytical treatment of Comptonization are essentially based on the solution of Kompaneets 
equation which describes the interactions between radiation field and thermal electrons 
(Kompaneets 1956). Due to the mathematical complexity, however, previous analysis of 
Comptonization depended on simplifications such as the non-relativistic approximation and 
therefore the results were only applicable to a relatively small range of photon and electron 
energies (e.g. Sunyaev & Titarchuk 1980). In recent years, Titarchuk (1994) developed a modified 
analytical technique which took into account the relativistic effect and Klein-Nishina corrections, 
thereby extending the previous work to wider ranges of temperature and optical depth of the 
plasma clouds from which Comptonized photons emerge. 
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The analytical method, however, have several intrinsic limitations. First, all analytical models 
are based on solving certain types of radiation transfer equations (Kompaneets 1956), which in 
turn is based on the assumption that energy and position of the photons are continuous functions 
of time, i.e. these models assume diffusion of photons in the energy and position spaces. While the 
continuity of energy change is a good approximation for scatterings at low energy, it is obviously 
not valid for Compton scatterings at high photon energies or by relativistic electrons. Similarly, the 
continuity of photon position change is an approximation only valid for clouds of electron plasma 
with dimensions large compared to the scattering mean free path (i.e. diffusion approximation). 
But astronomical observations suggest that many of the sources where Comptonization is believed 
to take place have optical depths of the order of one Thomson scattering mean free path. 

Second, solutions of the radiative transfer equations are based on the separation of photon 
diffusions in energy and position spaces (Sunyaev & Titarchuk 1980, Titarchuk 1994 and Hua 
& Titarchuk 1995). The solutions can be presented in terms of simple analytical expressions 
only when initial source photons have energies much lower than the electron energy and follow 
a particular spatial distribution, namely, the first eigenfunction of the spatial operator of the 
diffusion equation. It was found (Hua & Titarchuk 1995) that for source photons at energies 
not far below the electron energy or for clouds with large optical depth, the emergent spectra 
are sensitive to both the spectral and spatial distributions of source photons and the results of 
analytical method must be expanded to the higher order terms. Consequently, the analytical 
models are applicable only to certain ranges of plasma temperature and optical depth where 
solutions are insensitive to source conditions. 

Third, the analytical methods are inadequate to treat the temporal behavior of Comptonized 
emissions. Hua & Titarchuk (1996) have shown that for relativistic plasma, photons gain 
energy significantly with each scattering and consequently the scattering mean free path changes 
significantly with each scattering. Besides, for plasma clouds with small optical depth, the 
scattering mean free path are mainly determined by the boundary condition instead of the 
scattering cross sections. As a result, analytical treatment (e.g. Payne 1980), is only applicable 
to the limited situation in which electron plasma has non-relativistic temperatures and optical 
depths much greater than Thomson mean free path. 

In addition to the above limitations, analytical approach is totally incapable of dealing with 
the Comptonization problems involving complicated geometries and inhomogeneity of electronic 
media, where scattering mean free path depends on scattering location and direction as well as 
photon energy. But observations seem to indicate that investigations of Comptonization in the 
media with non-uniform temperature and density are necessary. As was shown by Skibo et al. 
(1995) and Ling et al. (1997), the spectral hardening at high energies in the spectra of AGNs and 
black hole candidates may be resulting from the temperature gradient in the plasmas responsible 
for the emissions. Kazanas et al. (1997) and Hua et al. (1997) showed that the temporal behavior 
such as the hard X-ray phase lags observed from the accreting compact objects may be explained 
by the non-uniform electron density of the accreting gas clouds. 
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These situations are where analytical method fails. As an alternative, Monte Carlo simulation 
can be employed to give solutions. It is flexible in simulating various initial conditions of source 
photons, complicated geometries and density profiles of plasma clouds. It is capable of presenting 
the full spectra resulting from Comptonization rather than the asymptotic ones obtainable from 
analytical methods. The first attempt to use Monte Carlo method to solve Comptonization 
problem was by Pozdnyakov et al. (1983). In recent years, Stern et al. (1995) presented a 
large-particle Monte Carlo method for simulating Comptonization and other high-energy processes. 
Skibo et al. (1995) used a Monte Carlo simulation in the calculation of photon spectra of mildly 
relativistic thermal plasmas in pair balance. 

In this study, we develop an efficient Monte Carlo method which treats Comptonization 
problem in a fully relativistic way and can be implemented in a medium computer such as 
Sparc workstation or Pentium PC to yield results with satisfactory statistics in CPU time of the 
order of minutes to hours. The algorithms are described in detail to facilitate computer code 
implementation. In §2 we introduce an improved technique of simulating Compton scattering of 
photons on cold electrons. In §3, we describe the method for Compton scattering on hot electrons. 
In §4, we present the method dealing with scattering in multi-zone medium. In §5, we describe the 
simulation of Compton scatterings in media with non-uniform density profiles. 



The Monte Carlo method described here was developed over the past several years in the 
investigations of Compton scattering of 2.223 MeV gamma-ray line in solar flares (Hua 1986), 
Compton backscattering of 511 keV annihilation line in the sources 1E1740. 7-2942 (Lingenfelter 
& Hua 1991) and Nova Muscae (Hua & Lingenfelter 1993). 

The differential cross section of Compton scattering is given by the Klein-Nishina formula 



where cjt is Thomson cross section; e = 2E/m e c 2 ; E is the energy of incident photon; m e is the 
electron rest mass and c the speed of light. The energy of the scattered photon, E', relative to the 
initial photon energy E is given by the ratio 



where ip is the angle between incident and scattered photons. The energy distribution of the 
Compton-scattered photons is determined by the distribution with respect to r, which is 



2. Compton Scattering on Cold Electrons 




(1) 



r = — = l + -(l-cosV>), 



(2) 




for 1 < r < £ + 1, 



otherwise, 



(3) 
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where 

K(e) = 3^a( £ ) (4) 

is the normalization factor. 

Sampling the distribution given by Eq. (3) plays a central role in the Monte Carlo simulation 
of Compton scattering of photons by cold electrons. Furthermore, as will be seen below, Compton 
scatterings on hot electrons in our scheme will also be reduced to the simulation of Eq. (3). 
Therefore, the performance of the computer program for Monte Carlo simulation of Compton 
scatterings depends critically on the quality of the technique used for sampling this distribution 
because a run of the program typically involves millions of scatterings. Efforts were made to 
optimize the technique of sampling this distribution (e.g. Kahn, 1954). In our implementation, 
we adopted a variation of Kahn's technique first suggested by Pei (1979). The algorithm of the 
technique is 



1. Generate 3 random numbers £1,62 and £3 uniformly distributed on (0,1). 

2. If £1 < 27/(2e + 29), 

let r= (e + l)/(e£ 2 + 1). 
If £3 > {[(e + 2 - 2r)/e} 2 + l}/2, go to 1. 
Else accept r. 
Else 

let r = e& + 1. 

If & > 6.75(r- l) 2 /r 3 , go to 1. 
Else accept r. 

It is seen that this is essentially a combination of composition and rejection methods (see e.g. 
Ripley, 1987). This algorithm, like Kahn's, avoids the operations such as square root, logarithm 
or trigonometric functions, which involve time-consuming series expansion for computers. Its 
quality can also be measured to a large extent by the rejection rates, which are 0.38, 0.30, 0.23 
and 0.33 for e = 0.2,2, 10 and 20 respectively, as compared to 0.41, 0.37, 0.41 and 0.53 for Kahn's 
technique. The improvement is significant, especially for higher photon energies. 



3. Comptonization in Hot Isothermal Homogeneous Plasmas 

The Monte Carlo technique for photon Comptonization in a relativistic plasma was outlined 
by Pozdnyakov et al. (1983) and Gorecki & Wilczewski (1984). Our implementation of the 
simulation is somewhat different from these authors. It was developed on the bases of the 
technique for Compton scattering on cold electrons described in the last section. 

Suppose a photon is scattered off an electron which is moving in z-axis direction with a 
velocity v. The energies of the incident and the scattered photon are E and E' respectively. The 
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zenith angles of the incident and scattered photons measured from z-axis are 9 and 9' respectively. 
4> and (f)' are the azimuthal angles. The differential cross section for Compton scattering is given 
by (see e.g. Akhiezer & Berestetskii 1969) 



da 3<tt 1 X ( E' s2 



(5) 



dfi'dft 16vr 7 2 (1 - vfi) 2 \E J ' 

where jx = cos 9 and // = cos 9'; v is in units of the speed of light and 7 = (1 — w 2 ) -1 / 2 ; 

e e' 4 / e\ 4 / e\ 2 

2<E IE 1 

' _ - e ' = — i^ 1 " V); (7) 



m e c* m e c 
E' 1-vfi 



E I — v/j,' + (E/jm e c 2 ) (1 — cos ip) 



(8) 



and tp is the angle between incident and scattered photons cos?/' = W t ' + V0- ~ ^ 2 )(1 ~~ I- 1 ' 2 ) cos(4> — 



Integration over // and 0' leads to 
3ut 1 



4 e 



1 - - - M 1 + £ ) + \ + - - 



e £ 2 J K ' 2 e 2(1 +e) 2 



(9) 



It is seen that Eq. (9) is identical in form with Eq. (1). But the quantity e here is given by the 
relativistic expression in Eq. (7). In other words, it is dependent on the electron's energy and 
direction as well as photon's energy. 

A photon with energy E traveling in a plasma with an isotropic distribution of electrons 
having an energy distribution ^(7) will have an averaged cross section of Compton scattering 
(see e.g. Landau Sz Lifshits, 1976): 

a a (T e , E) = -J d 7 J dn(l - vfi)a(£)N e ( 7 ). (10) 
For a plasma in thermal equilibrium, ^(7) is the Maxwell distribution given by 

where = kT e /m e c 2 is the dimensionless temperature of the plasma; k is the Boltzmann constant 
and K2 is the modified Bessel function of 2nd order. The <J a (T e , E) values in the form of a data 
matrix, obtained by the 2-dimensional integration in Eq. (10) for a properly spaced array of T e 
and E, can be read by or incorporated into the computer codes. Values of a a (T e ,E) for several 
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temperatures are numerically calculated and plotted in Figure 1. The dashed curve in the figure is 
the cross section at T e = 0, given by the Klein-Nishina formula in Eq. (1). It can be seen that for 
energetic photons scattering off the high temperature electrons, the cross section can be smaller by 
a factor of 2 or more than off the cold electrons. In other words, hot plasmas are more transparent 
than cold ones for photons. This has important effect on the energy spectra emerging from such 
plasmas, which Titarchuk (1994) took into account in his modification of the previous analytical 
results. Its effect on the temporal behavior of X-ray and gamma-ray emission from these plasma 
is even more significant and was discussed in Hua & Titarchuk (1996). 
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Fig. 1. — Maxwellian averaged Compton scattering cross section for various plasma temperatures, 
obtained from numerical integration in Equation (10). Also plotted is the maximum effective cross 
section as a function of photon energy H${E). 



With a a (T e ,E) obtained by numerical integration in Eq. (10), we can use the Monte Carlo 
method to select the free path between two successive scatterings for a photon with energy E. 



f 

Jo 



n e a a ds = -ln£, 



(12) 



where t is the free path to be sampled; n e is the electron density and £ is a uniform random 
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number on (0,1). The integration is taken along photon's path length s. In this section we are 
only concerned with the isothermal plasmas at temperature T e and with uniform density n e and 
leave the discussion about inhomogeneous media to the next two sections. Under this assumption, 
i can be sampled simply by 



n e a a (T e ,E)- 



(13) 



At the location of scattering, an electron is selected to scatter the photon. Its energy factor 
7 and direction [i = cos 9 with respect to the photon direction are selected according to the 
distribution 

(14) 



/ e ( 7 ,/x) oc (1 - ^l-^)a(e)N e ( 7 ), 



while its azimuthal angle (j) around the photon direction is selected uniformly on (0,27r). The 
distribution in Eq. (14) is rather complicated because e depends on 7 and /i as given in Eq. (7). 
On the other hand, for a thermal plasma, N e {^f) is given by Eq. (11) and independent of [i. In our 
implementation of the distribution Eq. (14), we use the following algorithm. 

1. Generate 2 random numbers £1 and £2 uniformly distributed on (0,1). 

2. if e < 0.01, 

If H > -e& I n 6, goto 1. 
Else let v = V-381nfi, 7 = - v 2 . 

Else if 6 < 0.25, 

let 7 = 1- 1.56 In ^1. 



If 6#i > Tv/Ci^ 2 -!), go to 1. 
Else 

let 7 = 1 - 361n£i. 

If 6#2 >7^i 2 v / 7^T, goto 1. 

3. Generate \x uniform on (-1,1) and £3 uniform on (0,1). 

4. Calculate e and then a(e) from 7 and {i according to Eqs. (7) and (9). 

5. if > (1 - ^i-rW, go to 1. 

Else accept 7 and /(/. 



Here 



Fi = aVa 2 - 1 



cxp 



a = 29 + 2b cos 



1 _! ^ 169 2 -1 

3 C ° S 



H 2 = aVa 2 - 1 



cxp 




a = + 2b cos 



1 _ x ^ 49 2 -l ' 
3 C ° S ^6^ 



1/3 + 49 2 ; 



(15) 
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1/3 + 9 2 ; 



(16) 
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and H3 is the maximum of the so called "effective cross section" a e s = (1 — n^/l — 7~ 2 )o"(e). 

Several points should be made in the above algorithm. Steps 1 and 2 sample 7 using rejection 
method in terms of the Maxwellian distribution iV e (7), which is independent of the photon energy 
and direction. For low plasma temperature (G < 0.01) electron velocity v are sampled according 
to the non-relativistic Maxwellian distribution. For high temperatures, the separated sampling 
(B < 0.25 and > 0.25) is in order to reduce the rejection rates. It should be emphasized that 
although the expressions of H\ and H2 are complicated, these quantities depend on alone and 
therefore need to be calculated once only. They can be calculated outside the scattering loop 
as long as plasma temperature remains unchanged. The 7 values so sampled, together with the 
isotropically sampled n, represent electrons in the hot plasma at the given temperature. They 
are subject to another rejection test in the subsequent steps in order to yield the right joint 
distribution given in Eq. (14), which represents the electrons that actually scatter the photon. 

temperature remains unchanged. The 7 values so obtained are subject to another rejection 
test in the subsequent steps together with the isotropically sampled \i in order to yield the right 
joint distribution given in Eq. (14). 

The quantity #3 is not expressible analytically. It depends on incident photon energy E only 
and can be determined by maximizing the effective cross section with respect to 7 and n for any 
given E using numerical methods such as given in Press et al. (1992). In the following, we describe 
an alternative to the above 2-dimensional maximization methods. We examine the derivative of 
o" e ff(7,/x) with respect to ji 



Therefore, 0^(7, /i) is a monotonously decreasing function of fi, that is, for given 7, <7 e ff(7,/x) 
reaches its maximum at \i = — 1. Physically, this means that head-on collision between the photon 
and electron always has the maximum probability. Thus, in order to determine the maximum of 
<r e ff as a function of 7 and /j,, one only needs to maximize the one dimensional function <r e ff(7, — 1). 
The maximum of <7 e ff, or H3, as a function of E determined in this way is plotted in Figure 1 
as the dash-dotted curve. It is seen that for high photon energies, the maximum effective cross 
section approaches the Klein-Nishina cross section while at low energies it approaches twice the 
Thomson cross section. The #3 values for an array of properly spaced E values can be tabulated 
and incorporated into the computer codes. 

With the selected electron energy and direction represented by 7, fi and </> uniform on (0, 2ir), 
we proceed to determine the energy and direction of the scattered photon. In order to do so, 




(17) 



where 




(18) 




(19) 
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we simulate Compton scattering in the frame where the electron before scattering is at rest 
rather than sampling the multivariate distribution of E', jjl 1 and <fi' from Eq. (5). The Lorentz 
transformation of the photon momentum between this reference frame and the lab frame is given 
by 

P' = P - Pbv - (7 - l)p • v]v, (20) 

where p and p' are photon momentum vectors before and after the transformation; p and v are 
unit vectors of the photon momentum and electron velocity respectively. 

In the electron rest frame, we utilize the Monte Carlo method described in §2. The resulting 
momentum of the scattered photon is then transformed back to the lab frame using the same Eq. 
(20) with a reversed v. The energy and direction of the scattered photon obtained in this way 
automatically satisfy the energy conservation relationship given in Eq. (8). 

As a crucial test we ran the program in which low frequency photons were allowed to scatter 
in an infinite plasma at a given temperature for a sufficiently long time. It was expected that the 
photon energy should approach the Wien distribution at the given plasma temperatures. One 
example of such evolution, the photon energy distribution recorded at varies times in a plasma of 
kT e = 200 keV are shown in Figure 2. It does approach the Wien form. 



4. Comptonization in Multi-Zone Media 

If Comptonization takes place in a medium which is divided into several zones each with 
different electron temperatures and density distributions, one has to take into consideration the 
boundaries between these zones in addition to the scattering free paths and the boundary of the 
entire medium. 

In general, suppose a photon, after initiation or scatterings, is located in the medium at 
(xo,yo,zo) with a direction (u±, 0J2, ^3). The next position where the photon will scatter, if there 
were no boundaries, is given by 

yi = y + £oj 2 (21) 
z\ = z + £lu 3 , 

where t is sampled according to Eq. (13), in which n e and T e should be understood as the electron 
density and temperature in the present zone. With the existence of boundaries, (x±, yi,z±) could 
be in the neighboring zone or outside of the medium. In this case, one has to calculate the 
distances Sj from (xo, Vo,zq) to various boundaries Bi,(i = 1, N), where ./V is the number of 
boundaries surrounding the zone under consideration. Si can be obtained by solving the equations 
describing the ith boundary 

Bi(x,y,z) = 0, i = l,...,N (22) 
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Fig. 2. — The evolution of photon energy spectrum from a blackbody at 0.511 x 1CT 2 eV towards 
equilibrium with a plasma at kT e = 200 keV. The seven spectra (solid curves) are "snapshots" at 
times t = 1,3,6,10,18,30,70 Thomson mean free time. Also plotted (dashed curve) is the Wien 
spectrum at temperature 200 keV. 



If i is smaller than any of si, ...sn so obtained, the photon will remain in the same zone and 
scatter at the location (xi,yi,zi) on electrons at local temperature T e . But if Sj is the minimum 
among t and s±, ...sat, the photon will hit the boundary Bj. In this case one has to replace the 
photon on the boundary at (x, y, z) determined by Eq. (23) with i = j. With the new position on 
the boundary as (xo, yo, zq), one can begin another round of free path sampling with n e and T e of 
the zone the photon is entering but keeping the photon energy and direction unchanged. 

In the study of Gamma-ray spectra of Cygnus X-l (Ling et al. 1997), we developed a model 
where photons scatter in a two-layered spherical plasma consisting of a high-temperature core 
and a cooler corona. The model was first proposed by Skibo and Dermer (1995) to interpret the 



where 




(23) 



X-ray spectral hardening at high energies observed in AGNs. The boundary of the inner core is a 
sphere of radius Ri while the boundaries of the outer shell are two spheres with radii Ri and R 
respectively. For a photon in the core, the equation for the distance si to its boundary is 



s\ + 2(r • Cj) Si - (R 2 - rl) = 0, 



(24) 



where ro = (xo,yo,zo) is the position vector of the photon; Co = (loi,lo2,cos)- Similarly, the 
equations for a photon in the outer shell are 



Thus we have the following algorithm: 
If r < Ri, 

Let S = (r ■ Co) 2 + (R 2 - rl) and si = y/S — (r ■ Co). 
If £ < si, scatter at ri = ro + £Co. 
Else reach boundary atri=ro + siu). 
Else if r < R , 

Let 5 = (r ■ Co) 2 + (R 2 - r 2 ). 
If 5 > and (r ■ Co) < 0, 

Let si = —s/S — (ro • Co). 

If I < s±, scatter at ri = ro + ICj. 

Else reach boundary atri = ro + sia). 
Else 

Let 5 = (r • Co) 2 + (R 2 - rl) and s 2 = V6 - (r • Co). 
If £ < S2, scatter at ri = ro + £Co. 
Else escape. 

Whenever the photon crosses the inner boundary, the plasma density and temperature should be 
switched while the photon energy and direction kept unchanged. 

In figure 3, we present the result of such a calculation (solid curve) together with the 
observational data (Ling et al. 1997) it was intended to fit. The data was from the blackhole 
candidate Cygnus X-l observed by the detector BATSE on board satellite Compton Gamma- Ray 
Observatory. The fitting spectrum was obtained from a calculation with the two-layer model 
described above, where temperature is kT e = 230 keV for the inner core and 50 keV for the outer 
shell. The two zones are assumed to have the same electron density and the inner core has a 
radius 0.36 in units of Thomson mean free path, while the outer-shell radius is 1.3. The initial 
photons have a blackbody temperature of 0.5 keV and injected into the medium from outside. For 
comparison the best fit one can achieve by a single-zone plasma model is also presented (dashed 
curve). The model consists of a plasma sphere of radius 1.35 at kT e = 85 keV. The reduced \ 2 



sj + 2(r • Co) Sl - (R 2 
s 2 + 2(r • Co)s 2 - {R 2 




(25) 
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value is 2.6 for the single-temperature model as compared to 1.0 for the double-layer model. It is 
seen that by adding a hot central core to the Comptonization medium, the fit to the high-energy 
part of the observed spectrum is significantly improved. 




lcr 1 10° io 1 io a io 3 iO 3 

E 7 (keV) 

Fig. 3. — The energy spectra resulting from the double-layer Comptonization media (solid curve) 
and singe-temperature sphere (dashed curve. Both spectra are intended to fit the observational 
data from the blackhole candidate Cygnus X-l (Ling el al, 1997). 



5. Comptonization in Isothermal Media with Non-Uniform Density 

The media we considered so far are uniform, at least regionally, in density. It was found 
necessary to investigate the Comptonization in the media with non-uniform density profiles 
(Kazanas et al, 1997). In this section, we present the treatment of two spherically symmetrical 
configuration commonly found in astrophysical environment, one with electron density varying 
as p' 1 and the other as /?~ 3 / 2 , where p is the distance from the sphere center. The latter case 
represents the density profile of a gas free- falling onto a central accreting object under gravitational 
force (e.g. Narayan & Yi 1994), while the former represents that of an accreting gas with viscosity 
due to the interaction between the gas and the outgoing photons (Kazanas et al, 1997). 
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With density n e varying along the photon's path length s, the integration in Eq. (12) should 
be written as 

I = / n e (s)a a ds, (26) 
Jo 



where the dependence of n e on s is given by 

mpo 



n e (s) = < 



Vq + s 2 + 2sr$v 

3/2 

k {rl + s 2 + 2sr z/) 3 / 4 



for p 1 profile, 
for p~ 3 / 2 profile, 



(27) 



where po is the radius of the sphere within which the density profiles break down; no is the 
electron density at this radius; v = (ro ■ w)/ro; ro is the photon's position vector originated from 
the sphere center and Co its travel direction. 



£-r u + J£ 2 + r 2 l + 2£r u 



Substitute n e (s) in Eq. (27) into Eq. (26) and we obtain the integration for p 1 profile 
/ = n p a a In 

r (l - v) 

Eq. (12) then becomes / = — ln£. Solving this equation for £, we obtain 

(1 + v)tf + 2ur] - (1 - v) 



(28) 



£ = r - 



2 V 



(29) 



where n = exp(— In £/ no po& a )- Once a uniform random £ is selected on (0, 1), £ is determined by 
Eq. (29). 



For p 3 / 2 density profile, the counterpart of Eq. (28) is 



(30) 



where F(ip,k) is the Legendre elliptic integral of the 1st kind; sin$ = y/l — u 2 ; ipo and (pg are 
given by 

fcos^o = (l + ng)- 1 ^ 
\con<pt = (l + ^ 2 )-V 4 , 



and 



(31) 



u 



cos $ 
sint? 

£ + r cos 
ro sin $ 



(32) 



-14- 



Substituting the integration into Eq. (12), we obtain 



^ / In ^ / In In £ / rn sin # . . 

F = F w,, -=) q —J -5- , (33) 

V2 V2 ™0P0^a V 2 / 9 o 

where the right-hand side is a function of known variables. Call it /(£, ro, ■#). Solve Eq. (33) and 
we obtain 

cosy*/ = cn(/,-^=), (34) 

where cn(/, k) is the Jacobian elliptic function, which is the inverse of the elliptic integral F((pg, k). 
Computer routines for both elliptic integral and Jacobian elliptic function are available in many 
mathematical libraries (e.g. Press et al, 1992). Finally, i can be obtained from Eqs. (31) and (32) 



£ = r sintf^cn- 4 (/, -J=) - 1 - r cos 0. (35) 

Once £ is available, one can use the algorithms described in the previous section to determine if 
the photon scatters, escapes or hits the boundary. 

We used £ given in Eqs. (29) and (35) to study the Comptonization in a two-layer spherical 
model similar to that in the last section but with the outer layer having a or p~ 3 / 2 density 
profile. Specifically, we assume the density in the outer shell is given by Eq. (27) with po = Ri 
and the density of the inner core is constant n + . It is found that the energy spectrum of the 
X-rays emerging from such system is different from a uniform sphere with the same optical depth 
(Kazanas et al, 1997). More importantly, with the decreasing density profiles, the outer layer, or 
the "atmosphere" can extend to a distance much greater than the size a uniform system with the 
same optical depth can do, giving rise to the time variation properties on a much greater time 
scale. 

As an example, we show in Figure 4 two light curves, or the time-dependent fluxes, for X-ray 
photons escaping from two such core-atmosphere systems, one with p~ l and the other with p~ 3 / 2 
density profile for the atmospheres. For both density profiles, the temperature is 50 keV in the 
atmosphere as well as in the core; the total optical depth is 2 in terms of Thomson scattering 
and the radius of the inner cores is assumed to be 2 x 1CP 4 light seconds. The core density n + 
is slightly different from each other: 1.6 x 10 17 and 1.68 x 10 17 cm -3 for and p -3 / 2 profiles 
respectively. For the outer atmospheres, no in Eq. (27) are 0.4 x 10 17 cm -3 and 1.68 x 10 17 cm~ 3 
respectively. As a result the radii of the systems are 1.01 and 2.63 light seconds respectively. 

Photons of a blackbody spectrum at temperature 2 keV are injected at the center into the 
system. The Comptonized photons in the energy range 10 — 20 keV are collected in terms of their 
escape time, producing the light curves displayed in the figure. It is seen that these light curves are 
power-laws extending to the order of seconds followed by exponential cutoffs. The indices of the 
power-law are roughly 1 and 3/2 respectively, which was explained in Kazanas et al (1997). This 
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Fig. 4. — The light curves resulting from the core-atmosphere models. The atmospheres have p 
and p~ 3 / 2 density profiles respectively. 

temporal behavior is greatly different from the light curves from a uniform system, which decay 
exponentially from the very beginning of the emissions (Hua & Titarchuk, 1996). In addition, 
for a uniform system of the similar optical depth and an electron density of the order of 10 16 or 
10 17 cm -3 , the characteristic decay time of the light curves will be ~ 1 millisecond. The implication 
of the prolonged power-law light curves resulting from the extended atmosphere models for the 
interpretation of the recent X-ray observational data is discussed elsewhere (Kazanas et al. 1997, 
Hua et al. 1997). 

6. Summary 

We have shown that analytical method has intrinsic limitations in dealing with Comptonization 
problem and Monte Carlo simulation provides a useful alternative. We have introduced an efficient 
Monte Carlo method that can solve the Comptonization problem in a truly relativistic way. The 
method was further expanded to include the capabilities of dealing with Comptonization in the 
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media where electron density and temperature vary discontinuously from one region to the other 
and in the isothermal media where density varies continuously along photon paths. In addition 
to the examples given above for its application, the method was also used in the investigation 
of Compton scattering of gamma-ray photons in the accretion disks near black hole candidates 
(Lingenfelter & Hua, 1991) and in the Earth's atmosphere and the spacecraft material (Hua & 
Lingenfelter, 1993). 

The author would like to thank R. E. Lingenfelter and R. Ramaty for their long-term support 
and encouragement in the past decade during which the technique described here was developed. 
The author also wants to thank J. C. Ling, L. Titarchuk and D. Kazanas for valuable discussions 
and NAS/NRC for support during the period of this study. 
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